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ABSTRACT 

We compare the evolution of morphology in X-ray cluster data and in clusters 
obtained in a simulation of flat CDM, normalized to the observed cluster abundance. 
We find that the evolution rate in model clusters is significantly higher than in data 
clusters. The results is stricking, as our cluster data are all just for z £ 0.2, but is not 
contrary to expectations, as is known that the CDM spectrum is too flat, around the 
cluster scale, and therefore induces the presence of too many substructures. The cluster 
simulations were run using a TREESPH code and therefore include hydrodynamical 
effects. The test we performed, which turns out to be very sensitive, starts from the so- 
called power ratios introduced by Buote and Tsai, and makes use of the 2-dimensional 
generalization of the Kolmogorov-Smirnov test. The discrepancy between data and 
model is however directly visible through linear regression techniques and using simpler 
statistical tests. 

Key words: Galaxies: clusters: general - galaxies: evolution - X-ray: galaxies - 
cosmology: theory - cosmology: simulations. 



1 INTRODUCTION 

The principal aim of this work is to study the evolution of 
morphology in hydrodynamical simulations of galaxy clus- 
ters, for flat CDM cosmological models, normalized in or- 
der to reproduce the observed cluster abundance. We also 
compare its trend with observational X-ray data and find 
a substantial discrepancy on the rate of evolution of sub- 
structures, which seem to evolve faster in models than in 
data. 

This result is not surprising, as the slope of CDM spec- 
trum, around the cluster scale, is known to be too flat. 
E.g., the extra power parameter F = 7.13T0 -3 (os/o^s) 10 '' 3 
(""8,25 are mass variances on 8,25 /i _1 Mpc scales; h = 
-ff/100 kms _1 Mpc~ 1 ) is significantly lower in data than in 
model predictions: for APM galaxies, Peacock and Dodds 
(1995) found T = 0.23 ± 0.04; for the Abell/ACO sample, 
Borgani et al. (1997) found T in the interval 0.18-0.25; on 
the contrary CDM models predict V ~ 0.4 . 

What is however stricking is the sensitivity that the 
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tests we performed on hydrodynamical simulations, seem to 
have in appreciating the spectral slope. High T values im- 
ply that too many substructures intervene in the dynamical 
growth of fluctuations, on cluster scales. An excess evolution 
can be therefore expected. But available data cover only a 
narrow redshift interval (0 < z < 0.2) and simulation out- 
puts shall be used on the same interval, reproducing data 
distribution. The point is then that different rates seem al- 
ready visible on such redshifts, which concern a fairly recent 
cosmic evolution. 

The tests we use start from the a statistical tool intro- 
duced by Buote & Tsai (1995), who defined the so-called 
power ratios IT"*', which essentially derive from a multipole 
expansion of X-ray surface brightness. In sec. 2 more details 
on n'-" 1 '' definition will be reported. Much work on n^was 
already done, both in order to derive them from data and to 
compare them with simulations (Buote & Tsai 1995, Buote 
& Tsai 1996, Tsai & Buote 1996, Buote & Xu 1997). In a 
recent work, Buote e Xu (1997) gave some plots where the 
evolution of cluster morphology in data and CDM simula- 
tions are compared; similar plots are given also by us. 

Here, however, we present two essential improvements: 
(i) Instead of using pure CDM N-body simulations, we use 
the outputs of a TREESPH code, therefore including hydro- 
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dynamical effects. A large number of massive clusters are 
treated in this way and hydrodynamical effects are found to 
be important in shaping the distribution of the X-ray emit- 
ting baryon component. This will be shown in some plots 
where DM and baryon distributions are compared; but there 
is also a back effect of the different baryon distributions on 
the evolution of the potential well, which cannot be appre- 
ciated directly from the plots, (ii) Furthermore, to compare 
data and simulations, we use an advanced statistical discrim- 
inator (Peacock 1982, Fasano & Franceschini 1987), which 
allows us to find quantitative estimates of the probability 
that the same process gives rise to observational and simu- 
lated statistical outputs. The discrepancy we find is however 
directly visible in suitable evolutionary plots and can be also 
appreciated using a simpler statistical test. 

The plan of the paper is as follows. In Sec. 2 power ra- 
tios are defined and their evolution is discussed. In Sec. 3 
we give some details about the model and simulation pro- 
cedure followed to define clusters and study their evolution. 
Power ratios are computed in Sec. 4 where we also perform 
the comparison between data and simulations. Results are 
discussed in Sec. 5, where some general conclusions are also 
drawn. Appendix A includes a brief outline of the PF2 sta- 
tistical test, and some discussion of the reasons why its ex- 
tension to three dimensions is premature. 



POWER RATIOS: 
EVOLUTION 



DEFINITION AND 



In this section we shall report results of Tsai & Buote (1996) 
and Buote & Xu (1997), which are the starting point of our 
arguments. The steps to work out the power ratios from a 
model cluster are the following ones: (i) We assume that 
the X-ray emission power is oc p b (r) and project it on a 
random plane 7r r , yielding a surface brightness E(x, y) (apart 
of a proportionality factor), (ii) We find the projected mass 
center (or centroid) O for T,(x, y); in respect to O, the surface 
brightness can be expressed as E(P, ip); (iii) we use such E to 
generate a pseudo-potential $(R, <p), by solving the Poisson 
equation: 



V 2 $ = T,(R,ip) 



(1) 



(again, a constant factor, in front of the r.h.s. - e.g., the 
gravitational constant G, - would turn out to be irrelevant), 
(iv) We expand $ in multipoles; from the coefficients of such 
expansion, the power ratios Y[( m ) (R) are obtained. More in 
detail: 



n (m) (.R) = iog 10 (p m /Po) 

where 

Pm(R) = ^(a™ +A) ,Po = [ao ln(P/kpc)] 2 



and 



= f dxx m+1 [ 
J Jo 

= f dxx m+1 [ 
J Jo 



d<p[T,(xR,(p)R ] cos(my) 



difi[T,(xR, ip)R 2 ] sm(m(p) 



(2) 



(3) 



(4) 



(5) 



Henceforth, the constants which were unspecified cancel out, 
when the ratio of P m 's is considered. Moreover, owing to the 
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Figure 1. Redshift distribution of data clusters for R 
0.4,0.8, l^ft^Mpc 



definition of the centroid O, Pi vanishes and IT 1 ' should 
then be disregarded. As we restricted the analysis to m < 4, 
which account for substructureson scales not much smaller 
than R itself, there are only 3 significant ir m )(m = 2,3,4). 

Therefore, during its evolution, a cluster moves along a 
line (evoJutionary track) in the 3-dimensional space spanned 
by such n^ m ^ 's, whose 2-dimensional projections will be con- 
sidered below. Starting from a configuration away from the 
origin, corresponding to a large amount of internal structure, 
a cluster evolves towards isotropization and homogeneiza- 
tion. This motion, however, does not occur with a steady 
trend, as sudden bursts of structure appear when further 
matter lumps approach the cluster potential well, to be ab- 
sorbed by it. However, the evolutionary track eventually ap- 
proaches the origin, and this can be more easily appreciated 
by averaging over the contributions of several clusters. 

Actual data, of course, do not show the motion of a 
single cluster along the evolutionary track. Different clus- 
ters, however, lie at different redshifts and, in average, can 
be expected to describe a succession of evolutionary mo- 
ments. Power ratios for model clusters are to be set in the 
3-dimensional space spanned by n^'s, taking each model 
cluster at redshifts distributed as for data clusters. 

The data set is the same of Buote & Tsai (1996). 
Among X-ray cluster images, taken with ROSAT PSPC 
(Pfeffermann et al. 1987), and such that the PSPC cen- 
tral ring contains a cluster portion whose radius exceeds 
400/i~ 1 kpc, they took those contained in the HEASARC- 
legacy database and belonging to the Ebeling (1993) or 
Edge et al. (1990) samples. Out of the 59 objects selected 
in this way, they could estimate n^for 44 of them at 
R = 0.8/i _1 Mpc and for 27 of them at R = 1.2fe _1 Mpc. 
PSPC data give the X-ray surface brightness T.x(R, <p), 
which is to be used in the same way as the E's obtained from 
models, to work out n (m) for R = 0.4, 0.8 , 1.2/i _1 Mpc. 

In fig. 1 we give the redshift distribution of data clusters, 
for the three values of R. As we shall discuss further in the 
next section, we have 40 simulated clusters to deal with. We 
take each of them, only at a single z value, producing 50 
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z distributions similar to real data. In order to work out 
power ratios, we need to treat each cluster with reference to 
a given plane: 50 random planes are selected for each cluster 
and one of them is considered for each z sequence. 



3 THE SIMULATIONS 

The main contribution of this work is given by the simu- 
lations we run, to work out IT" 1 ', and that will be how- 
ever used for several other aims. We started from a box 
with a side L — 200 /i~ 1 Mpc, where we evolved a CDM 
model with h = 0.5 and fif, = 0.06. The initial spec- 
trum was normalized in order that the cumulative cluster 
number density n(> M), for M = 4.2 /i _1 10 14 M©, yields 
N c i — n(> M)L 3 = 32, according to data on cluster 
abundance (White et al. 1993, Biviano et al. 1993, Eke et 
al. 1996). This corresponds to <jg = 0.52. The simulation 
was run from z in = 6.3 to z — using a P3M code with 10 6 
particles (mass m p = 1.78 • 10 13 M Q ). The choice of such z in 
value is dictated by the opportunity of restarting from the 
same conditions when hydrodynamic simulations are then 
performed. The PM part of the code makes use of 256 3 cells 
and a smoothing radius of 520 kpc was taken for gravita- 
tional forces. 

At z — clusters were then identified using a FoF al- 
gorithm; we considered as friends, particles at a distance 
< 0.2A (A is the mean particle separation). The 40 most 
massive clusters were then located; the lightest one is made 
by 34 particles and has a mass M m ~ 6 • 10 14 M Q . 

For each of these clusters an hydrodynamic simulation 
was then performed, using a TREESPH code. The hydro- 
dynamical part of the code is based on the SPH method 
(Hernquist & Katz 1989). The gravitational forces are solved 
with a hierarchical tree algorithm and this part of the code 
is based on the public treecode of L. Hernquist (Hern- 
quist 1987). In order to run the hdyrodynamical simula- 
tion we located the cluster center at z — and found all 
particles within a radius r2oo, where the cluster density is 
~ 200 (p) (background density). Such particles were traced 
back to Zi n , in the original simulation cube, and a smaller 
cube of side L c , enclosing all of them, was then located 
at the cluster center. We found the approximate scaling 
L c ~ 45Mpc x (r 20 o/3.43Mpc). 

A high resolution lattice of 22 3 particles was set in 
the smaller cube and their position was then perturbed, 
using the same initial conditions of the cosmological sim- 
ulations, implemented by additional waves to sample the 
increased Nyquist frequency. This technique is similar to 
the one adopted by Katz & White (1993), Navarro, Frenk 
& White (1995). We used two grids, for baryon and CDM 
particles, whose masses mb/mcDM = fif>/(l — fif>) • 

The cube of side L c is then set inside a greater cube of 
side 2L C , with parallel sides and equal center. Such greater 
cube is initially filled by 22 3 CDM particles, whose spacing is 
double in respect of inner cube particles, and whose mass is 
8 mco Af(l + fift). Then particles falling in the smaller inner 
cube are eliminated. Heavier particle positions are also per- 
turbed, according to the initial conditions of the cosmolog- 
ical simulation. Finally, we use all particles inside a sphere, 
of radius L c and with the same center as cubes, for the high 
resolution TREESPH simulation. 
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Figure 2. Density projections of a simulated cluster at several 
redishift. 

Since we take a constant number of grid particles, their 
masses vary according to L c ; for the most massive cluster 
we met (M 20 o ^ 10 15 M Q ), it was m gas = 3.4 1O 1O M . The 
softening parameter for gravitational forces (e) was set to 
e g = 80 kpc for gas particles and scales as oc m 1 '' 3 . We use 
a tolerance parameter 8 = 0.7, without quadrupole correc- 
tions. Viscosity was set as in Hernquist & Katz (1989) taking 
their parameters a = 1 and (3 = 2. 

The minimum timestep was 3 l(fyr for gas particles 
and twice as much for CDM. In average each cluster required 
~ 6 hours of CPU time to be evolved from Zi n to z = on a 
RISK 6000 elaborator. Outputs from TREESPH simulations 
were preserved at various redshifts. Among the smaller ones, 
we used z t =0.15, 0.10, 0.049 and (i = 1,..,4). 

In fig. 2 we give a typical example of projected pb out- 
puts, at a set of redshifts < 1.27. Already from such figure, 
it is clear that CDM clusters are rich of substructures and 
characterized by a fast substructure evolution. In fig. 3 we 
show CDM and baryon distributions at z = for a typical 
cluster. Such figure shows that the baryon distribution is sig- 
nificantly smoothed out in respect to CDM. This obviously 
reflect on n (m) 's, although we shall report no quantitative 
details on this point here. 



4 CLUSTER EVOLUTION FROM POWER 
RATIOS BEHAVIOUR 

Clusters are expected to evolve by accreting lumps of 
baryons and CDM onto a deep potential well. As time 
elapses, nearby lumps are gradually absorbed and cluster 
substructures tend to vanish. Such behaviour is clearly vis- 
ible also in fig. 2. In fig. 4, instead, we give 9 plots for 
the space containing the evolutionary tracks projected on 
3 planes and at three different _R's. Cross indicate the loca- 
tion of data clusters. Filled dots refer to a single selection 
of simulated clusters. The trend visible in this selection of 
simulated clusters is however a general feature: a linear re- 
gression gives place to the dotted straight line reported in 
the plots. Although real data do not show such a marked 
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Table 1. Angular coefficients and standard deviations for linear 
fits. Indeces d s refers to data and simulations. In parenthesis we 
give also simulation errors from 20 bootstrap resamplings. 
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Table 2. Kolmogorov-Smirnov test on angular coefficients for 
the sequence yielding fig. 4 and averaged on 50 redshifts sets. 
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Figure 3. Baryon and dark matter density projections of a sim- 
ulated cluster 



trend, a linear regression was performed also on them. The 
continuous lines show the results of such linear fits to data. 

Fig. 4 shows that model lines are steeper than data lines 
in all cases. This behaviour is true for all model sets we 
explored and indicates a faster evolution of model clusters 
in respect to the real ones. In Table 1, we report the values 
of the angular coefficients m of the two straight lines of each 
plot, together with 1 o errors Am, obtained from the linear 
regression. 

Another way of testing the stability of m amounts to 
performing a bootstrap-resampling of the simulation points. 
In fact, the 2-dimensional point distributions are rather ir- 
regular and a bootstrap approach tests the reliability of the 
linearity assumption, together with the actual m values. 

Bootstrap resampling was performed selecting clusters 
with the correct z distribution, but allowing data repeti- 
tion. All resamplings confirmed that model fits are steeper 
than data. Furthermore, bootstrap errors are systematically 
smaller than errors from linear regressions. Average errors 
from 20 bootstrap resamplings are reported in parenthesis 
in Table 1, aside of those from linear fits. 
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In Table 2 we report the the results of a Kolmogorov- 
Smirnov test applied to the line slopes of the single case 
considered in Fig. 4, and the values obtained by averaging 
on the set of 50 realizations. The figures are the probabil- 
ities that the different slopes we find in data and models 
can be due to the same process, taking into account the 
number of data points in use. This test reproduces the vi- 
sual feeling that the two straight lines, coming from data 
and simulations, have a systematically different slope. (Us- 
ing bootstrap errors as Am, besides of being scarcely signif- 
icant, would lead to even smaller probabilities.) This test, 
however, is somehow misleading and actual probabilities are 
substantially greater than so. The extra input that we should 
eliminate, to obtain the probability that the observational 
data set is fit by artificial data sets, is the requirement of 
linear behaviour. This can be done by applying the gener- 
alization of the Kolmogorov-Smirnov test to 2-dimensional 
data sets, proposed by Peacock (1982) and widely tested by 
Fasano & Franceschini (1989). We shall refer to this estimate 
as PF2 test and more details on it will be given in Appendix 
A. In Table 3 we give the values of the single case of Fig. 
4 and the averaged values of the PF2-probabilities that the 
observational data set and each artificial data set are real- 
izations of the same process, for each IT m pair and each R 
value. Such average PF2-probabilities are not so small as 
those reported in Table 2 and essentially correspond to a 
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Table 3. PF2— test on U m for the sequence yielding fig. 4 and 
averaged on 50 redshifts sets. 



n, - (R = D.4 tr 1 Mpc) 



n B - (R = 0.4 tr 1 Mpc) 




n 3 - (R = 0.4 tr 1 Mpc) 




II, - (R = 8 IT 1 Mpc) 





n, - (R = i.s tr 1 upo) 




n 3 - (R = i.a tr 1 Mpc) 



Figure 4. Power ratios and evolutionary tracks for data and 
simulated clusters at three aperture radii. Crosses and continuous 
lines refer to the observational sample, filled dots and dashed lines 
refer to one of the artificial sequences of 40 clusters, with similar 
z distribution. Plots obtained from all 50 cluster sequences are 
strictly analogous. 
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~ 2cr discrepancy. However, even this sophisticated test is 
unable to work out a final value which takes into account all 
the above partial results. There can be however little doubt 
left that such a sequence of values, where a faster evolution 
in model in respect to data is systematic, could hardly be 
due to observational sample biases or to disregarding pecu- 
liar physical effects in the simulations, and is more likely to 
indicate that the model is not adept to reproduce real data. 



5 CONCLUSIONS 

It has been known for several years that COBE normalized 
CDM yields a cluster abundance exceeding observation by 
a factor ~ 20. In order to obtain a fair number of clusters, 
the quadrupole data has to be lowered down to ~ 9/iK, 
at more than 3a's from the observational value. However, 
CDM models are still currently in use as reference models; 
in fact, 0CDM models, mixed models or A models, used 
to obtain fair fits of observational data on various scales, 
involve more parameters. 

According to our results, however, even though CDM is 
normalized to yield a fair number of clusters, we cannot ex- 
pect their morphologies to be adequately approached by the 
model. This output is related to the slope of the transfered 
spectrum, as already outlined in Sec. 1, and as is expressed 
by the values of the extra power parameter F. It is however 
important and fairly unexpected that the effects of such dif- 
ferent slope are already visible in the average evolutionary 
tracks, for a fairly narrow redshift range (0 < 2 < 0.2) as 
the one we had in observational data. 

In a forthcoming paper, based on work in progress, 
we shall extend this analysis to other cosmological models, 
whose transfered spectrum has a F value consistent with ob- 
servations. This is important, in order to test how sensitive 
the power ratios FL m are to the F value and to check whether 
this test depends only on F or, e.g., also on the substance 
of the model considered. 
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APPENDIX A - The PF2 test 

As is known, the Kolmogorov-Smirnov (KS) test is used 
to evaluate the probability that 2 distributions on a single 
variable x are produced by the same process. If N events 
are located at Xi (i = 1, .., N), let Sn(x) give the fraction of 
events located at Xi < x. Clearly Sn is discontinuous in each 
Xi, where it increases by 1/N. Let two sets of A^i and N2 
events yield Sn ± and Sjv 2 ; the KS statistics is then defined 
as the maximum distance between them: 



D 



?{any x } I S Nl (x) ~ Sn 2 (x)\ 



(al) 



and it can be shown that the probability that the two data 
sets are produced by the same process can be estimated, by 
using D, as follows: Let be 



l/v 2 = 1/N! + 1/N 2 



(a2) 



and u = D(y + 0.12 + 0.11/^). Then such probability reads: 



KS(u) 



-1)' / exp(2r 2 u 2 ) 



(a3) 



Various numerical routines exist which allow to evaluate 
such function with any desired approximation (see, e.g., 
Press et al. 1992). 

While for one-dimensional distributions the meaning 
and the building of D are simple (it tells us the maximum 
difference of integrated distributions), in the case of two- 
dimensional distributions the events are not ordered and 
there is no immediate way of telling where such difference is 
maximum. A further problem is set by the degree of correla- 
tion amongst data. A limiting situation occurs if they locate 
along a straight line on a plane, and then the KS test would 
be suitable again. 



The problem was considered in detail by Peacock (1982) 
and Fasano & Franceschini (1987), whose approach (PF2) 
we shall briefly report here. If N data are located on the 
2-dimensional plane x, y and have average x,y, let us first 
define & — Xi X and r\i = yi—y. Their correlation coefficient 
is then: 



(o4) 



(notice that, for r — 1, the points would actually lie on a 
straight line) . Then, if we have two sets of Ni and N% data, 
with correlation coefficients rj and r^, defined according to 
eq. (a4), we define 



(o5) 



as square average correlation. 

The basic point, however, amounts to define a 4-value 
extension of the function Sn, for each set of N data: the 4 
values are the fractions of data in each of the four quadrants 
set by the point x, y where Sn is evaluated. The PF2 statis- 
tics D will be the maximum difference between Sn ± and Sn 2 
values, ranging both over points and quadrants. Using such 
D, the square average correlation r, and v, defined according 
to eq. (a2), we then define 



= v 2 D/[v + ^/l-f 2 (0.25i/ - 0.75)] 



(06) 



and, using eq. (a3), we obtain the required probability. 

It ought to be outlined that there are some restrictions 
to the significance of these probability estimators. In gen- 
eral, they become safer as v increases. However, in the one- 
dimensional case, the estimate is good for v £ 10-20 (corre- 
sponding to Ni ~ N2 ~ 8). In the two-dimensional case, the 
restriction is already much stronger, as it must be v XL 100- 
400 (corresponding to Ni ~ N2 ~ 25), and the further con- 
dition that the probability found is ^ 20 % is to be added. 

Our sets of data and simulation points are therefore 
just adequate to meet the above conditions. This is impor- 
tant also in view of a step forward that could be considered: 
treating directly the 3-dimensional distributions in the pa- 
rameter space spanned by n'" 1 ' (m = 2,3,4), instead of 
their two-dimensional projections. An algorithm extending 
the above test to 3-dimensions is a straightforward gener- 
alization; the range of validity of such algorithm, instead, 
should be accurately tested. A large portion of the work of 
Fasano and Franceschini (1987) amounted to perform such 
checks, for the two-dimensional case. However, the trend 
found beween 1 and 2 dimensions, seems to indicate that, in 
three dimensions, we need 2 sets of data containing ~ 10 2 
points. As far as simulations are concerned, this requires 
just doubling our simulation set. Data clusters, instead, are 
signicantly below such limit. 

This paper has been produced using the Royal Astronomical 
Society/Blackwell Science F>TpX style file. 
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